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Abstract 

Alexandrov's Theorem states that every metric with the global topology and local geometry required 
of a convex polyhedron is in fact the intrinsic metric of a unique convex polyhedron. Recent work by 
Bobenko and Izmestiev describes a differential equation whose solution leads to the polyhedron corre- 
sponding to a given metric. We describe an algorithm based on this differential equation to compute 
the polyhedron to arbitrary precision given the metric, and prove a pseudopolynomial bound on its run- 
ning time. Along the way, we develop pseudopolynomial algorithms for computing shortest paths and 
weighted Delaunay triangulations on a polyhedral surface, even when the surface edges are not shortest 
paths. 

1 Introduction 

Alexandrov's celebrated theorem [Ale42| IAle05j characterizes the metrics of convex polyhedra. More pre- 
cisely, a convex polyhedron in Euclidean 3-space, viewed as a two-dimensional surface, induces an intrinsic 
metric on the (surface) points of the polyhedron: the distance between two such points is the length of the 
shortest path connecting them, restricted to lie along the polyhedron. We can divorce the intrinsic metric 
from the extrinsic embedding in 3-space, and Alexandrov's Theorem will tell us whether such an abstract 
metric could have come from a convex polyhedron. 

The intrinsic metric of a convex polyhedron has three obvious properties. First, the metric is polyhedral: 
every point but finitely many exceptions (the vertices) looks flat in the metric, meaning that it has a 
neighborhood isometric to a flat disk. Second, the metric is (locally) convex: every circle of radius r has 
circumference at most 2irr. Finally, treated as a topological space, the metric is homeomorphic to a 2-sphere. 

Alexandrov's Theorem [Ale42, AlcOSj says that these three necessary conditions are also sufficient: every 
convex polyhedral metric M homeomorphic to a sphere can be isometrically embedded as a convex polyhe- 
dron, meaning that its induced intrinsic metric is exactly M. Furthermore, the convex polyhedron is unique 
up to isometry of 3-space (an extension of Cauchy's Rigidity Theorem |Caul3[ [SR34] V Thus the extrinsic 
geometry can be reconstructed purely from the intrinsic geometry. 

Unfortunately, Alexandrov's proof is not constructive, suggesting an algorithmic problem: given a convex 
polyhedral metric homeomorphic to a sphere, find an isometric embedding as a convex polyhedron. More 
precisely, the polyhedral metric can be specified by a complex of triangles with specified edge lengths and 
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adjacency between triangles. It is easy to check that the given metric satisfies the three Alexandrov conditions 
(one even follows from our input representation). The goal is to find (approximate) coordinates for the 
vertices that (approximately) satisfy the edge-length constraints and convexity. 

One motivation for this problem is the problem of folding a given polygon of paper into precisely the 
surface of a convex polyhedron. There are efficient algorithms to find one or all gluings of a given polygon's 
boundary to itself so that the resulting metric satisfies Alexandrov's conditions DDLO02, L096J. These 
algorithms produce the desired convex polyhedral metrics homeomorphic to spheres, knowing from Alexan- 
drov's Theorem that they correspond to actual 3D polyhedra that can be folded from the polygon of paper. 
But without an algorithm for Alexandrov's Theorem, we do not know how to compute these polyhedra. 

Sabitov [Sab02, Sab96b, Sab96a, Sab98, IDO07] showed how to enumerate all the isometric mappings 
of a polyhedral metric as a polyhedron (not necessarily convex), which immediately leads to an algorithm 
for Alexandrov's Theorem. This algorithm has the distinction of being exact on a real RAM supporting 
polynomial root finding, which can also be implemented on a binary computer with a logarithmic dependence 
on the desired accuracy e. Unfortunately, the necessary polynomials have degree 2 e ( m ) for a polyhedron 
with m edges |FP05] . leading to an exponential running time. Without the convexity constraint, this 
exponentiality is unsurprising, because there can be exponentially many isometric mappings and hence 
exponential output size. But it is not known how to accelerate this algorithm in the convex case. 

The desire for a practical algorithm for Alexandrov's Theorem, "either a polynomial-time algorithm or a 
numerical approximation procedure", is posed as [DO07, Open Problem 23.22]. While the polynomial-time 
challenge remains open (and perhaps unlikely) , we come close in this paper by attaining a pseudopolynomial- 
time algorithm. Our work is based on recent progress on the other goal, a numerical approximation procedure. 

Namely, recent work by Bobcnko and Izmesticv [BI08 (building on work of Volkov and Podgornova 
|VP71| ) provides a new proof of Alexandrov's Theorem. Their proof describes a certain ordinary differential 
equation (ODE) and initial conditions whose solution contains sufficient information to construct the em- 
bedding by elementary geometry. The work in BI08] was accompanied by a computer implementation of 
the ODE [Scc06 , which empirically produces accurate approximations of embeddings of metrics on which it 
is tested. 

We describe an algorithm based on the Bobenko-Izmestiev ODE, and prove a pseudopolynomial bound 
on its running time. Specifically, call an embedding of a convex polyhedral metric M e-accurate if the metric 
is distorted by at most a factor 1 + e, and e-convex if each dihedral angle is at most 7r + e. Then we show 
the following theorem: 

Theorem 1.1. Given a convex polyhedral metric M homeomorphic to a sphere with n vertices, ratio S 
between the largest and smallest distance between vertices, and deject (discrete Gaussian curvature) be- 
tween E\ and 2ir — eg at each vertex, an e^-accurate eg-convex embedding of M can be found in time 
O (n 915 / 2 5 832 /(e 121 ef 45 el 17 )) where e = min(e 6 /nS, e 9 e\/nS & ). 

The exponents in the time bound of Theorem 1 1 . 1 1 are remarkably large. Thankfully, no evidence suggests 
that our algorithm actually takes as long to run as the bound allows. On the contrary, our analysis relies on 
bounding approximately a dozen geometric quantities, and to keep the analysis tractable we use the simplest 
bound whenever available. The algorithm's actual performance is governed by the actual values of these 
quantities, and therefore by whatever sharper bounds could be proved by a stingier analysis. 

To describe our approach, consider an embedding of the metric M as a convex polyhedron in R 3 , and 
choose an arbitrary origin O in the surface's interior. Then it is not hard to see that the n distances J-j = Ovi 
from the origin to the vertices Vi, together with M and the combinatorial data describing which polygons 
on M are faces of the polyhedron, suffice to reconstruct the embedding: the tetrahedron formed by O and 
each triangle is rigid in R 3 , and we have no choice in how to glue them to each other. In Lemma 13.21 below, 
we show that in fact the radii alone suffice to reconstruct the embedding, to do so efficiently, and to do so 
even with radii of finite precision. 

Therefore in order to compute the unique embedding of M that Alexandrov's Theorem guarantees exists, 
we compute a set of radii r — {ri}i and derive a triangulation T. The exact radii satisfy three conditions: 

1. the radii r determine nondegenerate tetrahedra from O to each face of T; 
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2. with these tetrahedra, the dihedral angles at each exterior edge total at most 7r; and 



3. with these tetrahedra, the dihedral angles about each radius sum to 2tt. 

In our computation, we begin with a set of large initial radii = R satisfying Conditions 1 and 2, and write 
k = for the differences by which Condition 3 fails about each radius. We then iteratively adjust the 

radii to bring n near zero and satisfy Condition 3 approximately, maintaining Conditions 1 and 2 throughout. 



The computation takes the following form. We describe the Jacobian I |^ ) , showing that it can be 



efficiently computed and that its inverse is pseudopolynomially bounded. We show further that the Hessian 



described by the Jacobian, with some step size only pseudopolynomially small, makes progress in reducing 
\k\. The step size can be chosen online by doubling and halving, so it follows that we can take steps of 
the appropriate size, pseudopolynomial in number, and obtain an r that zeroes k to the desired precision in 
pseudopolynomial total time. Theorem 11.11 follows . 

The construction of |BI08| is an ODE in the same n variables rj, with a similar starting point and with 
the derivative of r driven similarly by a desired path for k. Their proof differs in that it need only show 
existence, not a bound, for the Jacobian's inverse, in order to invoke the inverse function theorem. Similarly, 
while we must show a pseudopolynomial lower bound (Lemma [5T7]) on the altitudes of the tetrahedra during 
our computation, the prior work shows only that these altitudes remain positive. In general our computation 
requires that the known open conditions — this quantity is positive, that map is nondegenerate — be replaced 
by stronger compact conditions — this quantity is lower-bounded, that map's inverse is bounded. We model 
our proofs of these strengthenings on the proofs in |BI08| of the simpler open conditions, and we directly 
employ several other results from that paper where possible. 

One subroutine in our algorithm is of independent interest, which computes the weighted Delaunay 
triangulation on a polyhedral surface. Mitchell, Mount, and Papadimitriou MMP87 solve the related 
problem of computing the Voronoi diagram on a polyhedral surface. However, their algorithm assumes that 
the edges in the given triangulation of the surface are shortest paths in the metric, which does not hold in our 
setting. We show that their algorithm works for general triangulations, though the running time increases 
from polynomial to pseudopolynomial. Unfortunately, it seems difficult to dualize the weighted Voronoi 
diagram and obtain a weighted Delaunay triangulation, because in the weighted case, Voronoi cells can be 
empty and not contain their site. Fortunately, we show that a dual transform is possible in the unweighted 
case (even though Delaunay edges need not be shortest paths between their endpoints) , which lets us compute 
the unweighted Delaunay triangulation. Then we show that an incremental flip-based reweighting algorithm 
lets us put the weights back into the problem, while only requiring 0(n 2 ) flips and 0(n 2 lgn) time. 

The remainder of this paper supplies the details of the proof of Theorem 11.11 We give background in 
Section [2j and detail the main argument in Section [3l We bound the Jacobian in Section 0] and the Hessian 
in Section [5l Some lemmas are deferred to Section [6] for clarity. Finally, Section [7] describes how to compute 
weighted Delaunay triangulations on polyhedral surfaces. 

2 Background and Notation 

In this section we define our major geometric objects and give the basic facts about them. We also define some 
parameters describing our central object that we will need to keep bounded throughout the computation. 

2.1 Geometric notions 

Central to our argument are two dual classes of geometric structures introduced by Bobenko and Izmestiev 
in [BI08] under the names of "generalized convex polytope" and "generalized convex polyhedron" . Because in 
other usages the distinction between "polyhedron" and "polytope" is that a polyhedron is a three-dimensional 
polytope, and because both of these objects are three-dimensional, we will refer to these objects as "gener- 
alized convex polyhedra" and "generalized convex dual polyhedra" respectively to avoid confusion. 





is also pseudopolynomially bounded. It follows that a change in r in the direction of smaller k as 
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First, we define the objects that our main theorem is about. 

Definition 2.1. A metric M homeomorphic to the sphere is a polyhedral metric if each x £ M has an open 
neighborhood isometric either to a subset of R 2 or to a cone of angle less than 2ir with x mapped to the 
apex. The points falling into the latter case are called the vertices V(M) = {vi}i of M, and they must be 
finite in number by compactness. 

The defect Si at a vertex vi £ V(M) is the difference between 2tt and the total angle at the vertex, which 
is positive by the definition of a vertex. 

An embedding of M is a piecewise linear map / : M — > R . An embedding / is e- accurate if it distorts 
the metric M by at most 1 + e, and e-convex if f(M) is a polyhedron and each dihedral angle in f(M) is at 
most 7r + e. 

A perfect embedding of a polyhedral metric M is an isometry / : M — > K 3 such that f{M) is a convex 
polyhedron. Equivalently, an embedding is perfect if O-accurate and 0-convex. 

Alexandrov's Theorem is that every polyhedral metric has a unique perfect embedding, and our contribu- 
tion is a pseudopolynomial-time algorithm to construct e-accurate e-convex embeddings as approximations 
to this perfect embedding. 

Definition 2.2. In a tetrahedron ABCD, write ZCABD for the dihedral angle along edge AB. 

Definition 2.3. A triangulation of a polyhedral metric M is a decomposition into Euclidean triangles whose 
vertex set is V(M). Its vertices are denoted by V(T) = V(M), its edges by E(T), and its faces by F(T). 

A radius assignment on a polyhedral metric M is a map r : V{M) — > R + . For brevity we write for 
r(v t ). 

A generalized convex polyhedron is a gluing of metric tetrahedra with a common apex O. The generalized 
convex polyhedron P — [M, T, r) is determined by the polyhedral metric M and triangulation T giving its 
bases and the radius assignment r for the side lengths. 

Write Ki = 2tt — JZjk ZvjOviVk for the curvature about Ovi, and <pij = ZviOvj for the angle between 
vertices Vi , Vj seen from the apex. 

Our algorithm, following the construction in [BI08J, will choose a radius assignment for the M in question 
and iteratively adjust it until the associated generalized convex polyhedron P fits nearly isometrically in R 3 . 
The resulting radii will give an e-accurate e-convex embedding of M into R 3 . 

In the argument we will require several geometric objects related to generalized convex polyhedra. 

Definition 2.4. A Euclidean simplicial complex is a metric space on a simplicial complex where the metric 
restricted to each cell is Euclidean. 

A generalized convex polygon is a Euclidean simplicial 2-complex homeomorphic to a disk, where all 
triangles have a common vertex V, the total angle at V is no more than 2ir, and the total angle at each 
other vertex is no more than tt. 

Given a generalized convex polyhedron P = (ill, T, r), the corresponding generalized convex dual polyhe- 
dron D(P) is a certain Euclidean simplicial 3-complex. Let O be a vertex called the apex, A4 a vertex with 

OAi = hi = 1/r, for each i. 

For each edge ViVj £ E(T) bounding triangles ViVjVk and VjViVi, construct two simplices OAiAjnAijk, 
OAjAijkAju in D(P) as follows. Embed the two tetrahedra OviVjVk, OvjViVi in R 3 . For each i' £ {i,j, k, I}, 
place Aii along ray Ov^ at distance hy , and draw a perpendicular plane Py through the ray at A;< . Let 
Aijk, Aju be the intersection of the planes Pi, Pj, Pk and Pj, Pi, Pi respectively. By a standard computation in 
inversive geometry, A^k and Aju are on the respective perpendicular rays from O through ViVjVk and VjViVi, 
so OAiAjuAijk and OAjAij^Aju share the orientation of OviVjVk and OvjViVi because the two tetrahedra 
are together convex at edge ViVj by local convexity. 

Now identify the vertices Aijk, Ajki, Akij for each triangle ViVjVk £ F(T) to produce the Euclidean 
simplicial 3-complex D(P). Since the six simplices produced about each of these vertices Aijk are all defined 
by the same three planes Pi,Pj, Pk with the same relative configuration in R 3 , the total dihedral angle about 
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each OAijk is 2ir. On the other hand, the total dihedral angle about OAi is 2ir — Ki, and the face about At 
is a generalized convex polygon of defect m. Let 



_ hj -hj cos (f>ij 



be the altitude in this face from its apex Ai to side AijkAju. 

Definition 2.5. A spherical simplicial 2-complex is a metric space on a simplicial complex where each 2-cell 
is isometric to a spherical triangle. 

A singular spherical polygon (or triangle, quadrilateral, etc) is a spherical simplicial 2-complex homeo- 
morphic to a disk, where the total angle at each interior vertex is at most 2tt. A singular spherical polygon 
is convex if the total angle at each boundary vertex is at most tt. 

A singular spherical metric is a spherical simplicial 2-complex homeomorphic to a sphere, where the total 
angle at each vertex is at most 2w. 

The Jacobian bound in Section |4] makes use of certain multilinear forms described in [BI08 . 

Definition 2.6. The dual volume vol(/i) is the volume of the generalized convex dual polyhedron D(P), a 
cubic form in the dual altitudes h. 

The mixed volume vol(-, •, •) is the symmetric trilinear form that formally extends the cubic form vol(-): 

vol(a,&,c) = ^(vol(a + fe + c) - vol(a + b) - vol(6 + c) - vol(c + a) + vol(a) + vol(6) +vol(c)). 

The ith dual face area Ei {g(i)) is the area of the face around Ai in D(P) , a quadratic form in the altitudes 

g(i) = {hij}j within this face. 

The ith mixed area •) is the symmetric bilinear form that formally extends the quadratic form Ei(-): 

Ei(a, b) = ^{E t (a + b)- £,(o) - Ek(b)). 



Let TTi be the linear map 



/, s a hj hi cos (j>ij 



sin <pij 

so that iTi(h) = g(i). Then define 

Fi(a,b) = Ei(TTi(a),TTi(b)). 
so that Fi(h, h) = Ei(g(i), g(i)) is the area of face i. 

By elementary geometry vol(h, h 7 h) = | J2i hiFi(h 7 h), so that by a simple computation 

vol(a, b, c) = - ^2 a>iFi(b, c). 

i 

2.2 Weighted Delaunay triangulations 

The triangulations we require at each step of the computation are the weighted Delaunay triangulations used 
in the construction of BI08] . We give a simpler definition inspired by Definition 14 of Gli05 . 

Definition 2.7. In a polyhedral metric M with a radius assignment r, the weight of a vertex v is the square 
of its radius, so that h(v) — r(v) 2 . 

The power ir v (p) of a point p against a vertex v in a polyhedral metric M with weights w is pv 2 — w(v). 

The center C(viVjVk) of a triangle ViVjVk € T(M) when embedded in M 2 is the unique point p such that 
Kviip) = ^v (p) — n v k (p), which exists by the radical axis theorem from classical geometry. The quantity 
TT Vi (p) = Tr(viVjVk) is the power of the triangle. 
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A triangulation T of a polyhedral metric M with weights w is locally convex at edge ViVj with neighboring 
triangles v l VjV k ,v v i vi if 7T,,, (C^u,-^)) > n(viVjVk) and n Vk (C(vjViVi)) > ^(vjViVi) when ViVjVk,VjViVi are 
embedded together in R 2 . The two inequalities are equivalent by a lemma in classical geometry. The 
triangulation is strictly locally convex at UjUj if the inequalities hold strictly. 

A weighted Delaunay triangulation for vertex weights w or radius assignment r on a polyhedral metric 
M is a triangulation T that is locally convex at every edge. 

We describe in Section [7| an^ algorithm Polyhedral- Weighted-Delaunay to compute a weighted 
Delaunay triangulation in time 0(n 3 S , /eg). 

The radius assignment r and triangulation T admit a tetrahedron OviVjVk just if the power of ViVjVk is 
negative, and the squared altitude of O in this tetrahedron is —ir{viVjVk). The edge v^j is convex when 
the two neighboring tetrahedra are embedded in R 3 just if it is locally convex in the triangulation as in 
Definition 12.71 A weighted Delaunay triangulation with negative powers therefore gives a valid generalized 
convex polyhedron if the curvatures Ki are positive. For each new radius assignment r in the computation of 
Section|3]we therefore compute a weighted Delaunay triangulation and proceed with the resulting generalized 
convex polyhedron, in which Lemma 16.71 guarantees a positive altitude and the choices in the computation 
guarantee positive curvatures. 

2.3 Notation for bounds 

Definition 2.8. Let the following bounds be observed: 

1. n is the number of vertices on M. By Euler's formula, |-E(T)| and \F (T)| are both 0(n). 

2. £\ = mini &i is the minimum defect. 

3. 82 = minj(5i — ft,-) is the minimum defect-curvature gap. 

4. £3 = miriyg e(t) tfiij is the minimum angle between radii. 

5. £4 = maxi Ki is the maximum curvature. 

6. £5 = mm v . v . Vk£ F(T) ZviVjVk is the smallest angle in the triangulation. Observe that obtuse angles are 
also bounded: ZviVjVk < tt — ZvjViVk < it — £5. 

7. £6 is used for the desired accuracy in embedding M. 

8. £7 = (max; j L )/(mini — 1 is the extent to which the ratio among the Ki varies from that among the 
5i. We will keep £7 < £s/47r throughout the computation. 

9. £8 = miiii(27r — 5i) is the minimum angle around a vertex, the complement of the maximum defect. 

10. £9 is used for the desired approximation to convexity in embedding M. 

11. D is the diameter of M. 

12. L is the maximum length of any edge in the input triangulation. 

13. I is the shortest distance ViVj between vertices. 

14. S = max(£>, L) /£ is the maximum ratio of distances. 

15. do = miiipgM Op is the minimum height of the apex off of any point on M. 

16. di = mm v . Vj£ E(T) d(0,ViVj) is the minimum distance from the apex to any edge of T. 
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17. g?2 = mini r % is the minimum distance from the apex to any vertex of M. 

18. H = l/d ; the name is justified by hi — l/f*j < 1 / c?o ■ 

19. R = max.; r;, so 1/H < Ti < R for all i. 

20. T = HR is the maximum ratio of radii. 

Of these bounds, n, £i,£s, and S are fundamental to the given metric M or the form in which it is 
presented as input, and D, L, and I are dimensionful parameters of the same metric input. The values e$ 
and £g define the objective to be achieved, and our computation will drive £4 toward zero while maintaining 
£2 large and £7 small. In Section [5] wc bound the remaining parameters £3, £5, R, do, d±, and g?2 in terms of 
these. 

Definition 2.9. Let J denote the Jacobian ( f 51 ) . ., and H the Hessian ( a d K J ) . ., . 



3 Main Theorem 

In this section, we prove our main theorem using the results proved in the remaining sections. Recall 

Theorem 11.11 Given a polyhedral metric M with n vertices, ratio S ( the spread ) between the diam- 
eter and the smallest distance between vertices, and defect at least e\ and at most 2ir — £g at each ver- 
tex, an EQ-accurate Eg-convex embedding of M can be found in time O (n 915 / 2 S' 832 /(£ 121 £f 45 £g 17 )) where 
e = min(£ 6 /nS', EgE^/uS 5 ). 

The algorithm of Theorem 1 1 . II obtains an approximate embedding of the polyhedral metric M in I 3 . Its 
main subroutine is described by the following theorem: 

Theorem 3.1. Given a polyhedral metric M with n vertices, ratio S (the spread,) between the diameter and 
the smallest distance between vertices, and defect at least £\ and at most 2n — Eg at each vertex, a radius 
assignment r for M with maximum curvature at most £ can be found in time O (n 915 / 2 S' 832 /(£ 121 £f 45 £g 17 )) . 

Proof. Let a good assignment be a radius assignment r that satisfies two bounds: £7 < £s/47r so that 
Lemmas I6.5H6.7I apply and r therefore by the discussion in Section 12.21 produces a valid generalized convex 
polyhedron for M, and £2 = f2(£f Eg/n 2 ^ 2 ) on which our other bounds rely. By Lemma \Q. 11 there exists a 
good assignment r . We will iteratively adjust r° through a sequence r* of good assignments to arrive at an 
assignment r N with maximum curvature e^ < e as required. At each step we recompute T as a weighted 
Delaunay triangulation by algorithm Polyhedral- Weighted-Delaunay of Section [7] 

Given a good assignment r = r n , we will compute another good assignment r' = r n+1 with £4 — 
£4 = (£f 45 £l 21 £| 16 /(n 907 / 2 S' 831 )). It follows that from r° we can arrive at a satisfactory r N with N = 

O ( (n 907/2 5 831 )/(£ 121 £ 445 £ 61 6) )_ 

To do this, let J be the Jacobian (f^")^ and H the Hessian { d ^ K Q rk )^- fc , evaluated at r. The goodness 
conditions and the objective are all in terms of k, so we choose a desired new curvature vector n* in K-space 
and apply the inverse Jacobian to get a new radius assignment r' = r + J _1 (k* — k) in r-space. The actual 
new curvature vector k' differs from k* by an error at most \ |H||r' — r\ 2 < (||H| | J _1 | 2 ) \n* — k| 2 , quadratic 
in the desired change in curvatures with a coefficient 

by Theorems O and O and Lemmas |QRHllo"71 andlOl 

Therefore pick a step size p, and choose k* such that 

k* - m = -pKi - p ^Ki - Si min . (1) 
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Consider a hypothetical r* that gives the curvatures k* , and examine the conditions on £4, £2, £7 in turn. 

Both terms on the right-hand side of (HJ are nonpositive, so each m decreases by at least pKi. Therefore 
the maximum curvature £4 decreases by at least pe^. If any defect-curvature gap 5i — Ki is less than £i/2, 
then it increases by at least pKi > p(8i — £i/2) > p(e\/2); so the minimum defect-curvature gap £2 either 
increases by at least p£i/2 or is at least £i/2 already. Finally, the — pKi term decreases each Ki in the same 
ratio and therefore preserves £7, and the — p («j — Si mhij(Kj/Sj)) term decreases each ratio Ki/Si by p times 
the difference (Ki/Si — mhij^Kj /Sj)) and therefore reduces £7 by ^£7. Therefore k* would satisfy all three 
conditions with some room to spare. 

In particular, if we choose p to guarantee that each k\ differs from k* by at most ps^/2, at most pei/2, 
and at most p(£i/47r)(min i then this discussion shows that the step from r to r' will make at least half 
the ideal progress pe^ in £4 and keep £2, £7 within bounds. 



as required. Any smaller p will also produce a good assignment r' and decrease £4 by at least pe^/2 
proportionally. 

As a simplification, we need not compute p exactly according to Rather, we choose the step size 
p l at each step, trying first p t_1 (with p° an arbitrary constant) and computing the actual curvature error 
\k'-k*\. If the error exceeds its maximum acceptable value pe\ e^/lQir 2 then we halve p f and try step t again, 
and if it falls below half this value then we double p l for the next round. Since we double at most once per 
step and halve at most once per doubling plus a logarithmic number of times to reach an acceptable p, this 
doubling and halving costs only a constant factor. Even more important than the resulting simplification of 
the algorithm, this technique holds out the hope of actual performance exceeding the proven bounds. 

Now each of the N iterations of the computation go as follows. Compute a weighted Delaunay triangula- 
tion T* for r* in time 0(n 3 S/e s ) by algorithm Polyhedral- Weighted-Delaunay. Compute the Jacobian 
J* in time 0(n 2 ) using formulas (14, 15) in [BI08] . Choose a step size p*, possibly adjusting it, as discussed 
above. Finally, take the resulting r' as r t+1 and continue. The computation of k* to check p t runs in linear 
time, and that of r' in time 0(n") where uj < 3 is the time exponent of matrix multiplication. Each iteration 
therefore costs time 0(n 3 S/es), and the whole computation costs time 0(Nn 3 S/es) as claimed. □ 

Now with our radius assignment r for M and the resulting generalized convex polyhedron P with cur- 
vatures all near zero, it remains to approximately embed P and therefore M in IR 3 . To begin, we observe 
that this is easy to do given exact values for r and in a model with exact computation: after triangulating, 
P is made up of rigid tetrahedra and we embed one tetrahedron arbitrarily, then embed each neighboring 
tetrahedron in turn. 

In a realistic model, we compute only with bounded precision, and in any case Theorem 13.11 gives us only 
curvatures near zero, not equal to zero. Lemma 13.21 produces an embedding in this case, settling for less 
than exact isometry and exact convexity. 

Lemma 3.2. There is an algorithm that, given a radius assignment r for which the corresponding curvatures 
Ki are all less than e = O (mm.{e§/nS, £g£? /nS 6 )) for some constant factor, produces explicitly by vertex 
coordinates in time 0(n S/eg) an e^-accurate Eg-convex embedding of M. 



Since 



min^i > (maxKj)(min<5j/<5j)(l + £7) 1 > £4(£i/27r)/2 = £i£ 4 /47r 



and since 




(2) 
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Proof. As in the exact case, triangulate M, embed one tetrahedron arbitrarily, then embed its neighbors 
successively. Call the resulting configuration Q. The positive curvature will force gaps in Q between the 
tetrahedra, but since the curvature around each radius is less than e, the several copies of each vertex will 
be separated by at most neD. Now replace the several copies of each vertex by their centroid, so that the 
tetrahedra are distorted but leave no gaps. Call the resulting polyhedron P and its surface metric M' . 
The computation of a weighted Delaunay triangulation takes time 0(n 3 S/e%) by Algorithm Polyhedral- 
Weighted-Delaunay, and the remaining steps require time 0(n). We claim this embedding is £6-accurate 
and eg-convex. 

To show £6-accuracy, observe that since each copy of each vertex was moved by at most neD from Q to 
P, no edge of any triangle was stretched by more than a ratio neS, and the piecewise linear map between 
faces relates M' to M with distortion neS < e§ as required. 

Now we show eg-convexity. Consider two neighboring triangles ViVjVk,VjViVi in T; we will show the 
exterior dihedral angle is at least —eg. First, consider repeating the embedding with OviVjVk the original 
tetrahedron, so that OviVjVk,OvjViVi embed without gaps. This moves each vertex by at most neD, and 
makes the angle viViVjVk convex and the tetrahedron viViVjVk have positive signed volume. The volume of 
this tetrahedron in the P configuration is therefore at least —neD 3 , since the derivative of the volume in any 
vertex is the area of the opposite face, which is at always less than D 2 since the sides remain (1 + o(l))Z? in 
length. 

Therefore suppose the exterior angle Z.viViVjVk is negative. Then by Lemma l5.4l and Lemma l6.4| 

. 3[viViVjV k ][viVj\ (neD 3 )D 576nS 6 

2[v i v j vi][v j v i Vk\ (i 2 £5/4:) 2 e 2 

and since £2 > £i/2 at the end of the computation, ZviViVjVk > — e23QAnS & / e\ > — £g as claimed. □ 

We now have all the pieces to prove our main theorem. 

Proof of Theorem \l.l[ Let e = O (mm(es/nS, egef/nS 6 )) , and apply the algorithm of Theorem 13. II to obtain 
in time O (n 915 / 2 5' 832 /(£ 121 £| 45 £g 17 )) a radius assignment r for M with maximum curvature £4 < e. 

Now apply the algorithm of Lemma l3.2l to obtain in time 0(n 3 ) the desired embedding and complete the 
computation. □ 



4 Bounding the Jacobian 

Theorem 4.1. The Jacobian J = (f^)^- has inverse pseudopolynomially bounded by | JT 1 1 = O ( " 2 ^3^ 4 R) ■ 

Proof. Our argument parallels that of Corollary 2 in |BI08| . which concludes that the same Jacobian is 
nondegenerate. Theorem 4 of [BI08I shows that this Jacobian equals the Hessian of the volume of the 
dual D(P). The meat of the corollary's proof is in Theorem 5 of [BI08 , which begins by equating this 
Hessian to the bilinear form 6vol(/i, •, •) derived from the mixed volume we defined in Definition 12.61 So we 
have to bound the inverse of this bilinear form. 

To do this it suffices to show that the form vol(h,x, •) has norm at least ^(^/Jp,^) f° r au vectors x. 

Equivalently, suppose some x has |vol(/i, x, z)\ < \z\ for all z; we show \x\=0(f^fR). 

To do this we follow the proof in Theorem 5 of [BIOS] that the same form vol(h, x, ■) is nonzero for x 
nonzero. Throughout the argument we work in terms of the dual D(P). 

Recall that for each i, tt^x is defined as the vector {xij}j. It suffices to show that for all i 

,2 ^ ( n 3 T 3 2 n 2 T 2 . . \ 
\^ x \i = 0(- T —R 2 + -—R\x\ 1 ) 
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since then by Lemma 



,o An l2 ^ / n 4 T 3 2 n 3 T 2 . 

xU < — 7t max TT^a; U = O —5-5 — R H ^ — ° F l 

£3 * \ £26364 £26364 



and since Ixh < v^l^h and < a + implies X < -i/a + 6, IxU — O ( 71 ' 3 T Rj . Therefore fix an 
arbitrary i, let g = 7T;ft, and ?/ = Tr-jO;, and we proceed to bound \y\%- 

We break the space on which Ei acts into the 1-dimensional positive eigenspace of Ei and its [k — 1)- 
dimensional negative eigenspace, since by Lemma 3.4 of [BI08 the signature of Ei is (1, k — 1), where k is the 
number of neighbors of Wj. Write A+ for the positive eigenvalue and — E^ for the restriction to the negative 
eigenspace so that E~ is positive definite, and decompose g = <?-|_ + <?_ , y = y+ + y_ by projection into these 
subspaces. Then we have 

G k E t (g,g) = \ + g\ - ET(g_,g_) ± \ + g\ - G_ 

Ei(g,y) = X+g+y+ - Er( g _, y _) 
Y ± Ei(y, y) = X+y 2 + - Er( y _, y _) ± \ +y \ - y_ 

and our task is to obtain an upper bound on Y— = ET (y_ ,y~), which will translate through our bound on 
the eigenvalues of Ei away from zero into the desired bound on \y\. 

We begin by obtaining bounds on \Ei(g,y)\, G_, G, and Y. Since \z\ > \vo\{h,x,z)\ for all z and 
vol(/i, x, z) — J^j z jFj(h>, x), we have 



Further, det ( e^v'v)) < ® because Ei has signature (1,1) restricted to the (y,g) plane, so by 



Ei(g,y)\ = \Fi{h,x)\ < l. 

Ei{g,y) Ei(y,y)j 

Lemma 14.31 2 

Y = Ei(y,y) < — . 

62 

On the other hand — \x\i < XjFj(x, h) = y\ hjFj(x, x), so 

If, ,R 2 , \ f nT 



Y = Ei(y,y) = Fi(x,x) > -f ( (n - 1)— H + |x|i ) > - ( —R 2 + R\x\ 1 

hi V £2 / \ 6 2 

Now G = Ei(g, g) > 0, being the area of the face about Ai in D(P). We have \Ei\ = 0{n/e 3 ) by construction, 
so G,G_ < G + G- < \E,\\h\ 2 = 0{nH 2 /e 3 ) and similarly G = 0(nH 2 /e 3 ). On the other hand we have 
G = n(s 2 /R 2 ) by Lemma S31 

Now, observe that X +y+ g + is the geometric mean 

A+y+<?+ - V / (A + .g2)( A+2/ 2 ) = ^/( G + G _)(Y + Y_) 
and by Cauchy-Schwarz (y_ , g_) < G_ y_ , so that 



1 > Ei(y,g) > V(G + G_)(y + F_)- v/gIfI 



= y^I , S T= + y/G + G- } 



v/gTgT+v/gT v y/Y + Y- + y/YZ' 

If y > 0, it follows that 

y _ < 2^G + GZ _ Q ( n 2 T 2 R 2 



G V e 2 e 3 
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If Y < 0, then 



so 



+ G_ y/YZ 
y < 2 v /gTgT ; 2(-Y")(G + G_ 



G v G 
and because X 2 < a + bX implies X < y/a + b, 

y ^ 2^G + GZ y/2(-Y)(G + G-) 



G JG 



so that 



F_ = O (max (-Y)-^- J j = O (^R* + —mi 



In either case, using Y < R/e 2 and Lemma l4~4l we have 



n 3 T 3 „, n 2 T 2 



12/11 = vl + \y-\t < \E^\ ((Y + YJ) +YJ) = 0( -t—R 2 + T—mi 

and the theorem follows. □ 
Lemma 4.2. \x\ 2 < (4n/e§)maxj |7Tjx| 2 . 

Proof. Let i = argmax i and let Vj be a neighbor in T of Without loss of generality let Xi > 0. Then 

Xi — Xj COS tf>ij 1 — COS (f>ij . , . . . . . 

(7Tja;)i = > ^ , — - — - = x i t&n{(t> i j/2) > Xifcj 2 > \x oo£r 3 /2 

smcpij s:no, y 

and it follows that 

|7r,a;| > \iTix\oo > |sc|ooe3/2 > |ar|e3/2-\/n 
which proves the lemma. □ 

Lemma 4.3. F t {h,h) > e 2 /R 2 - 

Proof. The proof of Proposition 8 in BI08 shows that a certain singular spherical polygon has angular area 
Si — Ki, where the singular spherical polygon is obtained by stereographic projection of each simplex of P* 
onto a sphere of radius 1/r, tangent to it. The total area of the polygon is (5i — Ki)/r 2 at this radius, so 
because projection of a plane figure onto a tangent sphere only decreases area we have Fi{h, h) = area(P*) > 
(S i -K i )/rf>s 2 /R 2 . " " □ 

Lemma 4.4. The inverse of the form E$ is bounded by \E^ 1 \ = 0{n/e^). 

Proof. We follow the argument in Lemma 3.4 of |BI08j that the same form is nondegenerate. Let £j(y) be 
the length of the side between A4 and Aj in D(P) when the altitudes hij are given by y. Since Ei(y) = 
\ 12 j £j(y)Vj ^ follows that Ei(a, 6) = ~ J2j £j( a )bj- Therefore in order to bound the inverse of the form E,- L 
it suffices to bound the inverse of the linear map £. 

Consider a y such that \£(y)\ao < 1; wc will show j 2/ j 00 

= 0(n/e 4 ). Unfold the generalized polygon 
described by y into the plane, apex at the origin; the sides are of length lj (y) , so the first and last vertex are 
a distance at most < n from each other. But the sum of the angles is at least £4 short of 2ir, so this 

means all the vertices are within 0{nj£/C) of the origin; and the altitudes yj are no more than the distances 
from vertices to the origin, so they are also 0(71/64) as claimed. □ 

The proof of Theorem 14. II is complete. 
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5 Bounding the Hessian 

In order to control the error in each step of our computation, we need to keep the Jacobian J along the 
whole step close to the value it started at, on which the step was based. To do this we bound the Hessian 
H when the triangulation is fixed, and we show that the Jacobian does not change discontinuously when 
changing radii force a new triangulation. 

Each curvature Ki is of the form 2ir — J2j k-w v k £T / - v j^ ViVk ^ so m analyzing its derivatives we focus 
on the dihedral angles /.VjOviVk- When the tetrahedron OviVjVk is embedded in R 3 , the angle ZvjOviVk is 
determined by elementary geometry as a smooth function of the distances among 0,Vi,Vj,Vi-. For a given 
triangulation T this makes n a smooth function of r. Our first lemma shows that no error is introduced at 
the transitions where the triangulation T(r) changes. 



Lemma 5.1. The Jacobian J = (g:).- is continuous at the boundary between radii corresponding to one 



triangulation and to another. 

Proof. Let r be a radius assigment consistent with more than one triangulation, say with a flat face ViVjVkVi 
that can be triangulated by ViVk as ViVjVk, VkViVi or by VjVi as VjVkVi, viViVj. Since the Jacobian is continuous 
when either triangulation is fixed and r varies, it suffices to show that for neighboring radius assignments 
r + Ar, the curvatures k obtained with either triangulation differ by a magnitude 0(|Ar| 2 ), with any 
coefficient determined by the polyhedral metric or the radius assigment r. 

Embed the two tetrahedra OviVjVk, OvkViVi or OvjVkVi, OviViVj together in ]R 3 , with distances [Oi>j], etc., 
taken from r + Ar. Of the ten pairwise distances between the five points in this diagram, eight are determined 
by M or the radii and do not vary between the ViVk and VjVi diagrams. Since the angles ZvjOviVk, etc., 
are smooth functions of these ten distances, it suffices to show that the remaining two distances [UjUfe], [VjVr] 
differ between the diagrams by 0(|Ar| 2 ). Letting X denote the intersection of the geodesies ViVk, VjVi on the 
face ViVjVkVi, we have [viVk] in the ViVk diagram equal to [viX] + [Xvk], while in the VjVi diagram ViXvk form 
a triangle with the same lengths [wjX], and a shorter [vkVk]- The difference between [viVk] in the two 

diagrams is therefore the slack in the triangle inequality in this triangle ViXvk, which is bounded by 0(\ Ar| 2 ) 
since the vertices have moved a distance 0(|Ar|) from where r placed them with Vi,X,v^ collinear. □ 

It now remains to control the change in J as r changes within any particular triangulation, which we do 
by bounding the Hessian. 

Theorem 5.2. The Hessian H = (gf$^) ijk is bounded in norm by 0(n 5 ' 2 S 14 R 23 / {ela%d\ D 14 )) . 

Proof. It suffices to bound in absolute value each element g°.gp fc of the Hessian. Since Ki is 2ir minus the sum 

of the dihedral angles about radius r,, its derivatives decompose into sums of derivatives 9 f^f 28 where 
ViViVm € F(T). Since the geometry of each tetrahedron OviViv m is determined by its own side lengths, the 
only nonzero terms are where j, k £ {i, I, to}. 

It therefore suffices to bound the second partial derivatives of dihedral angle AB in a tetrahedron ABCD 
with respect to the lengths AB, AC, AD. By Lemma T5.6I below, these are degree-23 polynomials in the side 
lengths of ABCD, divided by [ABCD] 3 [ABC] i [ABD] i . Since 2[ABC],2[ABD] > {D/S)d u 6[ABCD] > 
d (D/S) 2 sine 5 , and each side is O(R), the second derivative is O (S u R 23 / {e^df D 14 )) . 

Now each element in the Hessian is the sum of at most n of these one-tetrahedron derivatives - g^'.g^" m , 

and the norm of the Hessian itself is at most n 3 / 2 times the greatest absolute value of any of its elements, 
so the theorem is proved. □ 

Definition 5.3. For the remainder of this section, ABCD is a tetrahedron and the dihedral angle ZCABD 
on AB. 




Lemma 5.4. 



sm 



e 



3 [ABCD] [AB] 
2 [ABC] [ABD] 
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Proof. First, translate C and D parallel to AB to make BCD perpendicular to AB, which has no effect 
on either side of the equation. Now [ABCD] = [BCD][AB}/3 while [ABC] = [BC][AB}/2 and [ABD] = 
[BD\[AB]/2, so our equation's right-hand side is ^c][bd} = sin ZCBD = sin#. □ 

Lemma 5.5. Each of the derivatives , g^j, g^ D is a degree- 1 polynomial in the side lengths of ABCD, 
divided by [ABCD][ABC] 2 [ABD] 2 . 

Proof. Write [ABC] 2 , [ABD] 2 as polynomials in the side lengths using Heron's formula. Write [ABCD] 2 as 
a polynomial in the side lengths as follows. We have 36[ABCD] 2 = det([AB, AC, AD]) 2 = det(M) where 
M = [AB, AC, A~b] T [A~B, AC, AD 1 } . The entries of M arc of the form u-v = ±(|u| 2 + \v\ 2 - \u- v\ 2 ), which 
are polynomials in the side lengths. With Lemma 15.41 this gives sin 2 9 as a rational function of the side 
lengths. 



Now ff = 9 a" 6 1 V 1 ~ sul2 ^ f° r anv variable x, so the square of this first derivative is a rational 
function. Computing it in SAGE [Ste08j or another computer algebra system finds that for each x € 
{AB, AC, AD}, this squared derivative has numerator the square of a degree-10 polynomial with denominator 
[ABCD] 2 [ABC] 4 [ABD] 4 . The lemma is proved. □ 

Lemma 5.6. Each of the six second partial derivatives of 8 in AB, AC, AD is a degree-23 polynomial in the 
side lengths of ABCD, divided by [ABCDf[ABC] i [ABD) i . 

Proof. By Lemma 15.51 each first partial derivative is a degree-10 polynomial divided by 
[ABCD] [ABC] 2 [ABD] 2 . Since [ABCD] 2 , [ABC] 2 , [ABD] 2 are polynomials of degree 6, 4, 4 respec- 
tively, their logarithmic derivatives have themselves in the denominator and polynomials of degree 5, 3, 
3 respectively in the numerator. The second partial derivatives therefore have an additional factor of 
[ABCD] 2 [ABC] 2 [ABD] 2 in the denominator and an additional degree of 13 in the numerator, proving 
the lemma. □ 



6 Intermediate Bounds 

In this section we bound miscellaneous parameters in the computation in terms of the fundamental param- 
eters n, S, £i,Ss and the computation-driving parameter £4. 

6.1 Initial conditions 

Lemma 6.1. Given a polyhedral metric space M , there exists a radius assignment r with curvature skew 
£7 < £s/47r, maximum radius R = OinD/exe^), and minimum defect- curvature gap £2 = £l(£ 2 £g/n 2 S 2 ). 

In the proof of Lemma 16.11 we require a lemma from singular spherical geometry. 

Lemma 6.2. Let C be a convex singular spherical n-gon with one interior vertex v of defect k and each 
boundary vertex Vi a distance a < W4 < f3 < n/2 from v. Then the perimeter per(C) is bounded by 

2,tt — k — 2n(n/2 - a) < per(C) < (2vr - n) sin /3. 

Proof. Embed C in the singular spherical polygon B that results from removing a wedge of angle k from a 
hemisphere. 

To derive the lower bound, let the nearest point on the equator to each Vi be Ui, so that UiVi < ir/2 — a. 
Then by the triangle inequality, 

per(C) = ViVj > u^j — v^i — UjVj > 2tt — n — 2n{ir/2 — a). 

ij ij 

For the upper bound, let D be the singular spherical surface obtained as the /3-disk about v in B. Then 
C can be obtained by cutting D in turn along the geodesic extension of each of the sides of C. Each of 
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these cuts, because it is a geodesic, is the shortest path with its winding number and is therefore shorter 
than the boundary it replaces, so the perimeter only decreases in this process. Therefore per(C) < per(D) = 
(2tt - k) sin /3. □ 



Proof of Lemma \6.1\ Let r have the same value R on all vertices. We show that for sufficiently large 
R = 0(nD/exEg) the assignment r is valid and satisfies the required bounds on £ 2 and £7. To do this it 
suffices to show that £ 2 < <5; — n% < £7£i for the desired £21 £ 7 an( i each i. 

For each vertex Vi , consider the singular spherical polygon C formed at Vi by the neighboring tetrahedra 
ViOvjVk- Polygon C has one interior vertex at ViO with defect Ki, its perimeter is Yljf. ^VjViVk = 2tt — Si, 
and each vertex ViVk is convex. The spherical distance from the center ViO to each vertex ViVk is ZOviVk = 
tt/2 — Q{viV k / R), which is at least p min = 7r/2 — Q(D/R) and at most p max = n/2 — Q(l/R). Now by 
Lemma RT21 above, we have 

2tt - Ki - 2n(ir/2 - p min ) < 2ir - 5{ < (2n - «,•) sin/9 max . 

The left-hand inequality implies 

S t - < 2ti(-k/2 - p min ) = 0{nD/R) 

so that Si — Ki < (£g/47r)£i if R = £l(nD/e\e%) for a sufficiently large constant factor. The right-hand 
inequality then implies 

Si -Ki> (2rr - ( 5 i ) 1 ~. SmPmax > £ 8 (1 - sin/w) = n(e 8 f/R 2 ) = n(eje 3 8 /n 2 S 2 ) 

sm Pmax 

so that the £2 bound holds. □ 

6.2 Two angle bounds 
Lemma 6.3. £ 3 > ldi/R 2 . 

Proof. £3 is the smallest angle <pij from the apex O between any two vertices ViVj. Now ViVj > £, and the 
altitude from O to v^j is at least d\. Therefore \ld\ < [OviVj] < \ smtfiijR 2 , so <pij > smcpij > ld\jR 2 . □ 

Lemma 6.4. £5 > E2/6S. 

Proof. Suppose that a surface triangle has an angle of e; we want to show e > £2/65. Let the largest angle 
of that triangle be tt — e'. By the law of sines, < S, so e > sine > sine'/S > d /3S since e' < 27r/3 
implies sine'/e' > 1/3. It therefore suffices to show that e 1 > £ 2 /2. 

Let the angle of size tt — e' be at vertex i. Embed all of the tetrahedrons around Ovi in space so that 
all the faces line up except for the one corresponding to an edge e adjacent to this angle of tt — e' . The two 
copies of e are separated by an angle of Ki. Letting / be the other side forming this large angle, the angle 
between one copy of e and the copy of / is tt — e'. Now the sum of all the angles around V{ is 2tt — Si, so 
apply the triangle inequality for angles twice to deduce 

£2 < 2tt - (2tt - Si) - ^ 

<2tt- ((tt - e') + Zfe') - k, 
= tt + e — Z/e — Ki 
< Ti" + e' - ((tt - e') - Ki) - Ki 
= 2e'. 

□ 
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6.3 Keeping away from the surface 

In this section we bound O away from the surface M. Recall that d 2 is the minimum distance from O to 
any vertex of M, d\ is the minimum distance to any edge of T, and d is the minimum distance from O to 
any point of M. 

Lemma 6.5. d 2 = f2(-D£i£4£§£ 8 /(nS' 4 )) . 

Proof. This is an effective version of Lemma 4.8 of |BI08j . on whose proof this one is based. 

Let i = argmin^ Ovi, so that Ovi = d 2 , and suppose that d 2 = 0(DeiS4e1es/(nS 4 )) with a small constant 
factor. We consider the singular spherical polygon C formed at the apex O by the tetrahedra about Ovi. 
First we show that C is concave or nearly concave at each of its vertices, so that it satisfies the hypothesis 
of Lemma \G. 91 Then we apply Lemma \6. 91 and use the fact that the ratios of the Kj are within £7 < £ 8 /4-7r 
of those of the Sj to get a contradiction. 

Consider a vertex of C, the ray Ovj. Let ViVjVk,VjViVi be the triangles in T adjacent to ViVj, and embed 
the two tetrahedra OviVjVk, OvjViVi in R 3 . The angle of C at Ovj is the dihedral angle VkOvjVi. 

By convexity, the dihedral angle VkViVjVi contains O, so if O is on the same side of plane VkVjVi as Vi 
is then the dihedral angle VkOvjVi does not contain Vi and is a reflex angle for C. Otherwise, the distance 
from O to this plane is at most Ovi = d 2 , and we will bound the magnitude of Z.V]~OvjVi. 

By Lemma I5T41 

. 3 [Ov k v 3 vi}[Ovj] 

sm Zv k OvjVi = -j- J -rj^ — J -\. 

2 [OvjV k ][OvjVi] 

Now [OvkVjVi] < d,2[vkVjVi]/3 = 0(g?2-D 2 ) and [Ovj] < [Ovi\ + [viVj] < D + d 2 - On the other hand [OvjVk] — 



(l/2)[Ovj][Ovk]sm ZvjOvk, and [Ouj], [Otik] > I — d 2 while Zv k ViVj < ZviV k O + Zv k Ovj + ZOvjV. 



ZvkOvj + 0(d 2 /D) so that ZvjOvk > £5 — 0{d 2 /D), so [OvjVk] = ^(^£5), and similarly [OvjVi] 
sin ZvkOvjVi — O (d 2 D 3 /(^ 4 £|)) = O (£i£4£s/n) , and the angle of C at Ovi is 



< 

Therefore 



ZvkOvjVi = O (£i£ 4 £ 8 /n) . 

On the other hand observe that per(C) = J2jk,ViW h eF(T) ZvjOvk < J2jk(^ v j v * v k + 0(d 2 /D)) 
+ 0{nd 2 /D). 

Now apply Lemma EU to deduce that 

per(C) ~ 



2tt 



Ki + 0(£i£4£s) > 1 



2tt 



zZ k i 



> 



2tt 



so that 



+ 0(e 4 e 8 ) > 



1 

2^ 



0(nd 2 / £l D) 



0(nd 2 /D] 



> (l + o(£ 8 )) 



= (l + o(e 8 )) 



1 



2tt V 

4tt - S. 



I mm — - 

3 6j 



2r: V 



3lti 
K 



> (1 + o(£ 8 ))(1 + £ 8 /2tt) fnun 



so that since Ki/5i = ^(£4), 



I mm — 

3 5,j 



> (l + 0(£ 8 ))- 1 (l + £ 8 /2 7 r) 



which for a small enough constant factor on d 2 and hence on the 0(£ 8 ) term makes £7 > £ 8 /(4tt), which is 
a contradiction. □ 
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Lemma 6.6. d a = n(e 2 4 e 2 5 dl/DS 2 ) = n(De\e\4el/ (n 2 S w )) . 

Proof. This is an effective version of Lemma 4.6 of |BI08j . on whose proof this one is based. 

Let O be distance d\ from edge ViVj, which neighbors faces ViVjVk,VjViVi € F(T). Consider the spherical 
quadrilateral D formed at O by the two tetrahedra OviVjVk, OvjViVi, and the singular spherical quadrilateral 
C formed by all the other tetrahedra. We will show the perimeter of C is nearly 2tt for small di and apply 
Lemma 16.81 to deduce a bound. This requires also upper and lower bounds on the side lengths of C and a 
lower bound on its exterior angles. 

In triangle OviVj, let the altitude from O have foot q; then Oq = d± while ViO,VjO > d%, so 
Z.VjVtq, ZviVjq = 0(d\/d2). Also, qvi,qvj > d 2 — d%, so q is at least distance {d^ — d\) sines from any 
of ViVk, vj~Vj,VjVi,viVi, and O is at least (d^ — d\) sine 5 — d\ = £l(d 2 £§) from each of these sides. 

Now ZviOvj = tt — 0(di/d<z) is the distance on the sphere between opposite vertices Ovi,Ovj of D, so 
by the triangle inequality the perimeter of D is at least 2ir — 0(di/d 2 ). Each side of C is at least £1{e§) and 
at most tt — Vlie^d^J D). 

In spherical quadrilateral D, the two opposite angles Zv^OviVi, ZviOvjVk are each within 0(d\/ e§d 2 ) of 
the convex ZvkViVjVi and therefore either reflex for D or else at least tt — 0{d\j £5^2). To bound the other 
two angles ZviOviVj, ZvjOvkVi, let the smaller of these be 9; then by Lemma [57 

tt — = G>(sm 0) — O [ - — = O 1 



x (e 5 d 2 D/S) 2 J V 4dl 

Now there are two cases. In one case, d\ = f^e^G?^/ DS 2 ). In the alternative, we find that each angle 
of D is at least tt — £4/2 and each angle of C at most tt — £4/2. In the latter case applying Lemma RT51 to C 
finds that 2tt - 0(di/d 2 ) = 2vr - fl(ele 5 d 2 /D) so that di = Q(ele 5 d%/D). 

In either case d\ = Cl(iain(e4eldl/DS 2 ,e 2 e 5 dl/D)) = ^l{e\e\d\j DS 2 ), and the bound on d 2 from 
Lemma [6751 finishes the proof. □ 



Lemma 6.7. 

a 1 £4 a 1 £4 \ \ _ n / £ i £ 4 £ 5 £ \ 



d = n min d^,^,^ 



Proof. This is an effective version of Lemma 4.5 of |BI08j . on whose proof this one is based. 

Let O be distance do from triangle ViVjVk € F(T). Consider the singular spherical polygon C cut out at 
O by all the tetrahedra other than OviVjVk- We show lower and upper bounds on the side lengths of C and 
lower bounds on its exterior angles, show the perimeter per(C) is near 2-7T for small c?o, and apply Lemma l6.8 
to derive a bound. 

The perimeter of C is the total angle about O on the faces of the tetrahedron OviVjVk, which is 2tt — 
O(do/di). Each side of C is at least ^(£5) and at most tt — VL{d\/D). 

Let 9 be the smallest dihedral angle of ZviOvjVk, ZvjOv^Vi, ZvkOviVj. Then by Lemma [5741 

«-e = o( s ^) = o(^l D 2 )=o fS2Ddo 



(diD/S) 2 J \ d\ 

Now there are two cases. If 9 < tt — £4/2, then it follows immediately that do = n(d 2 E4/(S 2 D)). 
Otherwise, 9 > tt — £4/2, so the interior angles of C are more than £4/2. Applying Lemma [6751 the perimeter 
per(C) is at most 2tt — D.{mh\{e\e^,e\di/ D)), so that d a = $7(min(di£4£^ 2 , <£j^ 2 £4-D -1 / 2 )). The bound on 
d\ from Lemma 16.61 finishes the proof. □ 



6.4 Lemmas in spherical geometry 

These lemmas about singular spherical polygons and metrics are used in Section [6.31 above. 



Lemma 6.8. Let a convex singular spherical polygon have all exterior angles at least 7 and all side lengths 
between c and 2tt — c. Then its perimeter is at most 2tt — 0(7 2 c). 
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Proof. This is an effective version of Lemma 5.4 on pages 45-46 of [BI08] , and we follow their proof. The 
proof in BI08 shows that the perimeter is in general bounded by the perimeter in the nonsingular case. 
In this case consider any edge AB of the polygon, and observe that since the polygon is contained in the 
triangle ABC with exterior angles 7 at A, B its perimeter is bounded by this triangle's perimeter. Since 
c < AB < 2ir — c, the bound follows by straightforward spherical geometry. □ 

Lemma 6.9. Let S be a singular spherical metric with vertices {i>i}i, and let C be the singular spherical 
polygon consisting of the triangles about some distinguished vertex Vq. Suppose C has k convex vertices, each 
with an interior angle at least tt — e for some e > and an exterior angle no more than tt. Then 

«o + 2e*>(l-^)5>. 

Proof. We reduce to Lemma 5.5 from [BJ08 by induction. If k = 0, so that all vertices of C have interior 
angle at least tt, then our statement is precisely theirs. 

Otherwise, let v^ be a vertex of C with interior angle tt — 9 G [it — s,tt). Draw the geodesic from v^ to 
Vq, and insert along this geodesic a pair of spherical triangles each with angle 6/2 at Vi and angle kq/2 at 
Vq, meeting at a common vertex v . The polygon C" and triangulation S' that result from adding these two 
triangles satisfy all the same conditions but with k — 1 convex vertices on C", so 



4 + 2e(k -!)>(!- Per(C ' 



2tt J V 

Now C and C have the same perimeter, k,' q < kq + 6 < hq + e, — — 9 > «j — e, and k'j — Kj for 
j £ {0,?}, so it follows that 

kq + 2ek >k' + {2k - l)e > f 1 - ^ ) Yl Ki 
as claimed. □ 



7 Delaunay Triangulation on a Polyhedral Surface 

In this section we present an algorithm Polyhedral- Weighted-Delaunay to compute a weighted De- 
launay triangulation on a polyhedral surface, as defined in Section 12.21 above, and we prove its correctness 
and efficiency Our algorithm consists of three main parts. First, in Section [7.1l we compute the unweighted 
Voronoi diagram on a polyhedral surface, using a generalization of the Mitchell Mount -Papadimitriou al- 
gorithm to allow the surface to have non-shortest-path edges. Second, in Section [7T21 we use this Voronoi 
diagram to compute an unweighted Delaunay triangulation, which is complicated by the fact that there are 
many possible paths between two vertices. Finally, in Section [7.31 we show how to modify this triangulation 
into a weighted Delaunay triangulation, by continuously reweighting the vertices and performing flips as 
necessary. 

7.1 Shortest paths on a non-shortest-path triangulation 

In this subsection we describe modifications to the analysis of the "continuous Dijkstra" algorithm of 
MMP87] and |Mou85j , which can be used to compute an (unweighted) Voronoi diagram on a polyhedral sur- 
face. Our modifications permit the algorithm to dispense with the assumption that the input triangulation 
consists of shortest paths, at the cost of a modest loss in efficiency when the assumption is not satisfied. 

We first sketch the main ideas of |Mou85j and |MMP87j . These papers describe an algorithm for com- 
puting shortest paths on the surface of a polyhedron from a number of sources. The shortest paths on 
the surface are represented by the shortest paths to each edge, approaching through each adjacent face. 
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Specifically, each edge, for each adjacent face, is partitioned into intervals on which the shortest paths to 
each point through that face originate at the same source and pass through the same sequence of edges and 
vertices. Then for each such interval the algorithm considers the last vertex in the sequence and records the 
location of that vertex in an edge unfolding along the sequence of edges. It is easy to see that once this 
representation is obtained, standard planar techniques suffice to efficiently obtain the shortest path to any 
point and the Voronoi diagram on the surface. 

The authors describe their algorithm for computing this representation as a "continuous Dijkstra" algo- 
rithm because it proceeds in a fashion analogous to Dijkstra's algorithm for shortest paths in a graph. The 
MMP algorithm proceeds by maintaining in a priority queue the shortest yet-known paths to a number of 
intervals. At each step it removes the nearest interval in the queue, identifies the shortest yet-known path 
as an actual shortest path to at least the nearest point in the interval, and propagates the paths through 
the interval to one or more opposite edges of the next face. The output of the algorithm is correct by an 
induction, and the runtime of the algorithm is governed by a bound on the number of intervals it must visit 
and propagate. 

Two modifications are required in order to extend the MMP algorithm to suit our purpose. First, the 
definition of a Voronoi diagram in |Mou85] must be slightly modified: it includes a point x in the Voronoi 
cell of source s if x is as close to s as to any other source. A better definition includes x only if it has a 
unique shortest path to s, shorter than any path to another source. It is straightforward that the Voronoi 
cells under the latter definition are simply connected. Because the algorithm concerns shortest paths to edge 
points only, and the Voronoi diagram is constructed by standard planar techniques after the MMP algorithm 
proper is complete, this modification requires no change to the MMP algorithm itself. 

Second, and more complex, both |Mou85j and MMP87J make an assumption that the polyhedral surface 
is embedded in K 3 , or at a minimum that a lesser property holds which is not satisfied by a general poly- 
hedral metric represented by a general triangulation. The earlier paper [Mou85l explicitly disclaims such 
an assumption, but at the outset of Section 3 it asserts that the restriction of a shortest path to a face is a 
single line segment, which is indeed easy to see if the triangulation derives from a Euclidean embedding or 
otherwise if each edge is a unique shortest path, but is not true in general. The use of this assertion is in the 
complexity analysis, so we repeat the complexity analysis after proving a relaxed version of the assertion, 
consisting of a bound on the number of line segments that may make up the restriction of a shortest path 
to a face. This re-analysis makes up the remainder of this section. 

Our argument requires an adjustment to one concept used throughout the analysis in |MMP87j . For an 
edge e and a face / bordered by e, the original analysis describes the paths associated with the edge-face pair 
(e, /) as "/-free", leading from the source to e while never passing through /. In our case where faces need 
not consist of shortest paths, the same algorithm will consider paths which are not always /-free, but they 
will be f -facing in the sense that they lead into / when extended through their endpoint at e. The analysis 
in [MMP 8 7] goes through unchanged with this broader definition, with the exception of the complexity 
analysis. We rehearse the latter with our modifications after supplying ourselves with a series of geometric 
lemmas. 

The main object of our geometric study in this section will be the piecewise- geodesic loop, or simply loop, 
a non-self-intersecting closed path on M consisting of finitely many geodesic segments. By the Jordan curve 
theorem, each loop partitions the rest of M into two sides, and we will sometimes identify one side as the 
"interior" and the other as the "exterior". By a closed side of a loop we mean one side together with the 
loop itself. 

Lemma 7.1. If a piecewise- geodesic loop of length L has at least two vertices on each closed side, then 
L = Q(£e & ). 

Proof. We shall repeatedly "cut" an angle of the loop, as follows. Draw a chord of the loop very close to 
the angle, so that it encloses a triangle with no vertices inside, and scale up by a homothety until the chord 
meets either a vertex or an angle; if it meets an angle, continue sliding the other endpoint of the chord until 
it meets another angle or the chord meets a vertex. The essential feature of the cutting process is that it 
always decreases the total length of the loop (by the triangle inequality) , so that every loop we work with is 
of length at most L. 
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Now, cither the loop passes through two vertices, or it passes through one and has at least one vertex on 
each side, or it has two on each side. Clearly if it passes through two vertices it has length at least 21 and 
we are done. We shall first reduce the case of at least two vertices on each side to the case of one vertex on 
the loop and at least one on each side. 

At least one side of the original loop has total defect at most 2ir, so identify one such side as the "inside" . 
We shall move the loop consistently toward the inside. At least some angles of the loop open toward the 
inside, unless the loop has no angles at all, in which case it is a closed geodesic, the surface is locally a 
cylinder, and we may slide the loop perpendicularly until we hit a vertex and thereby create an angle. Now 
we cut an angle that opens toward the inside, which either eliminates a vertex from the interior, making it 
an angle of the loop, or reduces the number of angles in the loop. Therefore repeating this finitely many 
times brings us to a loop which contains just one vertex v in its interior, passes through at least one other 
vertex u, and still has length at most L. 

Once we reach this stage, we proceed by cutting any angle which is not the vertex u. If any such 
cut hits a vertex w, then u and w are within distance L/2 and L > 2£. Otherwise, we cut until u is 
the only angle remaining in the loop, and because we encountered no vertices there is still only the one 
vertex v in the interior. Therefore the interior is metrically a cone, and by elementary geometry L > 
2dist(w,u)sin(7r- S v /2) > 2^sin(e 8 /2) = ft(£s s ). □ 

Call a piecewise-geodesic loop live if it meets the condition of Lemma 17.11 having at least two vertices 
on each side. 

Consider a digon xGyHx, where x and y are points and G and H two geodesies. By local geometry, the 
extensions of G and H through x are on the same side of the digon, as are the two extensions through y. 
If the extensions through x and through y are on opposite sides, we call xGyHx an inside-out digon, and 
otherwise all four extensions are on the same side and we call it a normal digon. 

We will repeatedly use the following basic fact: on a polyhedral metric a digon always encloses at least 
one vertex on each side. Further, because by Gauss-Bonnet a normal digon encloses a total defect less than 
2ir, its closed exterior contains defect greater than 2ir and therefore at least two vertices. 

Lemma 7.2. // non-self-crossing geodesies G and H in M form an inside-out digon xGyHx with only one 
vertex on one side, then G and H end on that side after crossing at most 0(1/ Sg) times. 

Proof. Call the side with one vertex the interior, arbitrarily, and by renaming let y be the vertex through 
which G and H enter the interior. Extend G into the interior until its next crossing z with H. 

If z lies on xHy, then xGzHx and zGyHz are normal digons whose interiors partition the interior of 
xGyHx. But each of these interiors must contain a vertex, making at least two vertices inside xGyHx, a 
contradiction. 

If yGzHy is a normal digon, then it must enclose a vertex in its interior, and the region interior to 
xGyHx and exterior to yGzHy has no vertices. Now G and H continue into this region; extending G, it 
cannot cross itself, but if it crosses any segment of H it divides this region into two pieces, of which one is 
a digon, which with no vertices is impossible. So G must end in this region, and similarly H must end, and 
there is only one crossing z in the interior of xGyHx. 

Otherwise yGzHy is an inside-out digon, dividing the interior of xGyHx into a quadrilateral region 
xGyH zGyHx and a region bounded only by yGzHy, which we call the interior. This interior is an entire 
side of yGzHy, so it must contain the vertex. 

We have one crossing for each successive inside-out digon that G and H form, and at most one for a 
normal digon at the end of the geodesies. Now consider the angles at which G and H cross at each successive 
crossing x,y,z,.... The difference between successive angles is the total external angle of the interior of 
the digon, which equals 2tt minus the total defect enclosed, which is at most 2tt — eg. Therefore successive 
crossing angles increase by at least e$, so because each angle is less than it we have at most 7r/e 8 inside-out 
digons before the paths end. □ 

The next lemma is the one employed in our extended complexity analysis. 
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Lemma 7.3. On the polyhedral metric M , a non-self-crossing geodesic G of length L > £ and a shortest 
path H can cross at most 0(L/£es) times. 

Proof. Consider the sequence along G of its crossings with H. We show that for each crossing y, either its 
two neighbor crossings x and z are separated by a length at least Q,(te%) of G, or one side of y along G has 
at most 0(1/ 'eg) crossings in total. The desired bound follows. 

Let three consecutive crossings with H along G be x, y, and z. Then the two digons xGyHx and yGzHy 
partition M into three regions. 

If either digon is inside-out, say xGyHx, then either it is live or it has only one vertex on one side. In 
the latter case, by Lemma I7T21 there are 0(1/ Eg) crossings toward that side before G and H end, and we are 
done. 

Otherwise both digons are normal. Therefore they each have at least two vertices in their exterior, so 
if either one has at least two vertices in its interior then it is live. Otherwise they each have exactly one 
interior vertex, so because M has at least four vertices there are two vertices in the remaining region and its 
boundary, either the quadrangle xGyHzGyHx or the digon xGzHx if xyz are out of order on H, is live. 

Now we have a live digon, or a live quadrangle whose edge set is a union of digons. By Lemma 17. II the 
live digon or quadrangle has length VL(ie%). Each digon consists of one segment from G and one from H, so 
since H is a shortest path at least half the total length of the digon must be in G. Therefore we have at 
least a length Q,(£es) in either xGy, yGz, or their union, so that the length of the segment xGz is at least 
n(feg) a s claimed. □ 

Now we proceed to the argument of [MMP87 , making the necessary extension to the complexity analysis 
of the continuous Dijkstra algorithm. 

Lemma 7.4. At the conclusion of the continuous Dijkstra algorithm, each edge-face pair (e, /) has at most 
0(n 2 S/e%) intervals in its interval list. 

Proof. We follow the proof of Lemma 7.1 in MMP87J, extending it where necessary. First, consider a single 
source v, and let I\, 1%, . . . , Ir be the interval list for the edge-face pair (e, /) that would be produced from 
the single source s. Let x% be a point interior to interval and let xq and xk+i be the endpoints of e. As 
in |MMP87] . the intervals are ordered so that / is on the left as we pass from Ij to Ij+i, and we should 
name xq and xk+i so that / is on the left as we walk from x toward Xr+i- Draw the shortest /-facing 
paths from s to each Xi, calling them Pi. Now between each Pi and Pj+i we have a region sPiXiexi+iPi+is, 
the region on the right as we traverse that cycle. Clearly there is a vertex in the interior of that region, 
or else Pi and P,+i would intercept the same sequence of edges so that X{ and x i+ i would belong to the 
same interval. But the whole cycle sPqXqbxr+iPr+xS can wind about any particular vertex v at most 
0(D/£e$) — 0(S/e$) times, because a shortest path from v to a point on the outside will cross e only that 
many times by Lem.ma l7.3l Therefore at most 0(S/ss) of the K — 1 such regions may enclose any particular 
vertex, and so K = 0(nS/e%). 

Now to complete the proof, we study the actual interval list produced by the algorithm when it runs with 
all our sources. Each interval I is a contiguous region on which the shortest paths lead from a particular 
source s through a particular sequence of edges. Suppose that two intervals /, I 1 derived from the same 
source and sequence of edges, so that some other interval J intervened which derived from another source 
t, or from t — s but with a different sequence of edges. Then we may place s and t into the plane of / by 
edge-unfolding along the respective sequences of edges, so that the points of I and I 1 are closer to s than 
to t but the points of J are closer to t than to s. But this is impossible because our distances are decisive. 
Therefore at most one interval in the actual interval list derives from a given source and sequence of edges, 
so the interval list contains at most one interval for each interval in the per-source interval lists, for a total 
of 0(n 2 S/e$) intervals. □ 

Lemma 7.5. The continuous Dijkstra algorithm runs in time 0(n 3 S / e%) . 

Proof. By Lemma I7T41 there are at most 0(n 3 S/es) intervals in the algorithm's entire final list. Each call to 
Propagate causes one final interval, and creates at most two new intervals. Therefore there are 0(n 3 S/ss) 
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intervals ever created. The ordered-set operations required on interval lists take time 0(\og(nS/es)), there 
are O(l) of these and O(l) other work in each round of the algorithm, one round per final interval, and so 
the total time is 0(n 3 (S/ss) log(nS/e$)) = 0(n 3 S/es) as claimed. □ 



7.2 Unweighted Delaunay triangulation on a polyhedral surface 

Next we give an efficient algorithm for computing an unweighted Delaunay triangulation on a polyhedral 
surface. An algorithm based on successively flipping edges which fail the local convexity condition was 
previously known to terminate, but is not believed to finish in polynomial time [ILTCOlj IRiv94) . Instead, 
here we use the unweighted Voronoi diagram computed from the previous section. Note that the continuous 
Dijkstra algorithm can be modified to compute weighted Voronoi (power) diagrams, but it seems difficult 
to transform such a diagram into a corresponding weighted Delaunay triangulation; thus we focus on the 
unweighted case for now. 

Algorithm Polyhedral-Delaunay. Begin by computing the unweighted Voronoi diagram using the 
generalized MMP algorithm described in Section 17.11 For each two Voronoi cells x and y that touch, and 
each segment e on their mutual boundary, take a point p on that segment. Because p is on the two Voronoi 
cells, the disk D p of radius px — py about p has no vertices in its interior, so it is isometric to a planar disk. 
Let e' be the segment xy through D p . It is straightforward to see that e' does not depend on the choice of p. 

From the data provided by the MMP algorithm at p, we may compute in constant time the length of the 
radii px and py of D p and the angles they make at p, x, and y, so by trigonometry we may also compute in 
constant time the length of e' and the angles it makes at x and y. Let T consist of the segments e! for each 
Voronoi boundary segment e, described by their lengths and their ordering about each vertex. 

If the Voronoi diagram has any points p at which more than three cells meet, then the disk D p about 
p has the source of each cell on its boundary and no vertices in the interior. The chords between adjacent 
vertices on the boundary of D p are already present in T as the segments derived from the Voronoi edges 
meeting at p, and we add to T further chords chosen to triangulate arbitrarily the polygon they form. □ 

Lemma 7.6. Algorithm Polyhedral-Delaunay computes an unweighted Delaunay triangulation in O(n) 
time plus the 0(n 3 S/e%) time to compute the unweighted Voronoi diagram. 

Proof. In the conversion from Voronoi to Delaunay, we spend O(l) for each Voronoi edge, and 0(d) work 
for each Voronoi vertex of degree d, for a total of 0(n) work. 

For correctness, we need to show that the edges of T form a triangulation, and that the triangulation is 
locally convex at each edge. 

First, we show that the edges of T do not intersect. Suppose that edges e' and /', which were respectively 
drawn through disk D p from A to C and through D q from B to D, intersect at X. Now D p contains no 
vertices in its interior, so /' must extend from X at least to the boundary of D p in each direction, and by 
the classical theorem on the power of a point, XB ■ XD > XA ■ XC. Similarly XA ■ XC > XB ■ XD, so the 
two powers are equal and D p = D q is a disk with A, B, C, and D all on its boundary. Consequently p = q 
is a point at which at least the four cells for A, B, C, and D meet, and e' and /' are diagonal chords chosen 
in the final step of the algorithm. But these chords are chosen as a triangulation, so they do not intersect. 

Now, we count the edges of T and apply Euler's formula to deduce that they form a triangulation. The 
Voronoi diagram has n faces; let it have mi edges and mo vertices, so that mi — mo = n — 2 by Euler's 
formula. Each Voronoi vertex of degree d contributes d — 3 edges in the last step of the algorithm, so a total 
of 2mi — 3m edges come from this step, for a total of 3mi — 3m = 3n — 6 edges in T. Since T has n 
vertices, by Euler's formula it must have 2 + (3n — 6) — n — 2n — 4 faces. Therefore the average degree of a 
face is 2 • (3n — 6)/(2n — 4) = 3. Since M is connected, every set of Voronoi cells must share edges with its 
complement, so T is connected, and therefore by geometry it can have no digons or empty loops. Therefore 
every face of T is a triangle. 

Finally, consider an edge e' of T separating triangles ACD and CAB. In the unweighted case, the local 
convexity condition reduces to a requirement that when ACD and CAB are developed together into the 
plane, D is not in the interior of the circumcircle of CAB. But in this development, B and D lie on opposite 
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sides of line AC, and the disk D p through which e' was drawn touches A and C and contains neither B nor 
D in its interior. Therefore the circumcircle of CAB has its center at least as far to the B side of AC as D p 
does, and contains a subset of those points to the D side of AC that D p does, so that it does not contain D 
in its interior and e' is locally convex. □ 

7.3 Weighted Delaunay triangulation on a polyhedral surface 

Algorithm Polyhedral- Weighted-Delaunay finds a weighted Delaunay triangulation for a polyhedral 
surface M with vertex weights w satisfying certain conditions. The weights w must obey a system of linear 
inequalities that guarantee that any edge AC which is not strictly locally convex lies in a convex quadrilateral 
consisting of two neighboring triangles ABC and CDA, and therefore can be "flipped" to substitute the 
contrary diagonal BD. Proposition 4 of |BI08j guarantees that every w arising from a generalized convex 
polyhedron, which includes every w we consider, satisfies this condition. 

The algorithm makes use of the following geometric observation. In polyhedral surface M with vertex 
weights w and triangulation T, let the height function on M be the unique function h : M — > R that within 
each triangle of T is quadratic with unit quadratic part along every line segment, and such that at each 
vertex v the height agrees with the weight, h(v) = w(v). We call this function the "height" because when 
the vertex weights arise from a radius assignment r , it coincides with the squared distance from the apex 
in the tetrahedra built from T and r. (In this setting, h is the function gr,r of [BI08j.) Then we have the 
following lemma, which is mentioned in [BI08 without proof. 

Lemma 7.7. If T is not locally convex at edge xz with neighboring triangles xyz and wxz, then replacing 
xyz and wxz with xyw and wyz strictly increases the height of the points in the interior of wxyz. 

Proof. Consider the points in the intersection U of triangles xyz and xyw. This is without loss of generality, 
because every point inside wxyz lies in at least one equivalent such intersection. Develop both triangles into 
the plane, and for clarity of discussion, orient the diagram so that edge xy is horizontal, and the alternate 
vertices w, z are above xy. (Because xz is not locally convex, wxyz is a convex quadrilateral so that w and 
z are necessarily on the same side of xy.) 

The height function on segment xy is determined in both triangulations as the unique quadratic function 
with unit quadratic part along xy and agreeing with the vertex weights at x and y. The height function on 
U is easily seen to be determined by its gradient at xy, and the gradient of h at any point in a triangle ViVjVk 
is twice the displacement vector from the center C(viVjVk). Both centers Cixyz) and C(xyw) lie on the 
radical axis axis(x,y) of x and y, the line on which ir x (p) = ir y (p), which is perpendicular to xy. Therefore, 
the height function at any point in U \ xy is greater in xyz, greater in xyw, or equal in both triangles just if 
C(xyz) is lower than C(xyw) along axis(x, y), is higher than C(xyw), or coincides with C(xyw) respectively. 

But if T with triangles xyz and wxz fails to be locally convex at xz, then ir w (C(xyz)) < ■K x (C(xyz)), so 
that C(xyz) falls on the w side of the radical axis axis(x, w). Because w is above x and bxXs{x,w) _L xw, 
the w side is the upper side, so that C(xyz) falls above the intersection of &xis(x,w) with axis(x, y), which 
is C{xyw). Therefore the height function is greater in triangle xyw than in triangle xyz, as required. □ 

Let XY be an edge of the triangulation, and let P be on XY. Then 

py.m ) + Pi., ( y) 

V ; XY 



Therefore in a pair of triangles WXY, WYZ forming a quadrilateral WXYZ with P = WY n XZ, if WY 
is in the triangulation we have 

v ' WY 
and if XZ is in the triangulation we have 

h{P) = P^^) + Px. w{Z )_ px pz 
X z 
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so that for any assignment of weights to W, X, Y, and Z we may determine by a comparison of these 
quantities whether triangles WXY and WYZ would flip to XYZ and XZW or vice versa. In particular, 
consider w(X), w(Y), and w(Z) as fixed. Then WXY and WYZ are preferred over XYZ and XZW or 
vice versa just if w(W) is greater than or less than 

^ . WY ( PZ.«X) + PXMZ) _ PW^Xl _ PX . PZ + PW . „.) , (3) 

Now, suppose we have a weighted Delaunay triangulation T for some vertex weights w, and consider a 
new vertex weighting w' which differs from w at only one vertex v, with w'(v) < w(v). The local convexity 
criterion is identical in w' and w for all pairs of adjacent triangles neither of which is incident to w. Further, 
by our discussion of equation |3J any edge in T which lies across a triangle from v remains locally convex 
under w' because the weight of v has only declined. 

The following algorithm Reweight makes use of T to compute a Delaunay triangulation T for w' . 

Algorithm Reweight. Let v have degree d. For each three consecutive neighbors x,y,z of v in T, define 

ty = tv.xyz- Compute the d values t y using equation [3] and store them in a max-priority queue Q. Now 
repeat the following steps until done: 

Remove the maximum element t y from Q. If t y < w'(v), we are done. Otherwise, replace edge vy with 
xz. Recompute t x and t z to reflect the modified triangulation, and update the priority queue. 

When the iteration is done, the resulting triangulation is T'. □ 

Lemma 7.8. Algorithm Reweight converts a Delaunay triangulation for w into a Delaunay triangulation 
for w' in time O(nlogn). 

Proof. There are d < n elements in the priority queue, each of which is inserted once and removed at most 
once, and there are two updates and one operation on T for each removal. This totals 0(d) operations which 
may each be done in time O(logd), for a total time 0(d\ogd) = 0(n log n). 

We prove correctness by induction. Let Reweight change m edges, let w(v) = to > t\ > ■ ■ ■ > t m+ i 
with ti the value removed from Q in step i, and let T = To, T%, . . . , T m — T be the successive triangulations 
considered. Write w t for the vertex weighting that differs from w only in setting the weight of v to a value t. 
We claim that for i = 0, . . . ,m, triangulation Tj is Delaunay for weightings w t with U >t> U+i. 

For i = 1, . . . ,m, triangulation T,-_i is Delaunay for w ti by hypothesis. Triangulation Tj differs only in 
one pair of edges, at which by construction it remains locally convex for w ti . Now by the discussion above, 
the only edges in Tj that may fail local convexity under any wt with t < U are those incident to v. But tj+i 
was constructed as the maximum weight for v at which the local convexity condition would reach equality 
for any of these edges. Therefore Tj remains Delaunay for weightings Wt with U > t > tj + i as claimed. 

In the base case of i = 0, triangulation To = T is Delaunay for w by precondition. Then To remains 
Delaunay for weightings Wt with to > t > t\ by the same argument as in the inductive case, and the induction 
is complete. 

Now setting i = m, we find that T = T m is Delaunay for weightings w t with t m > t > t m +\. But by 
the termination condition, t — w'(v) lies in this range. Consequently T is Delaunay for us w i/ v > = w' as 
required. □ 

With Polyhedral-Delaunay and Reweight as subroutines, it is now straightforward to compute a 
weighted Delaunay triangulation for any vertex weighting w on a polyhedral surface M. 

Algorithm Polyhedral- Weighted-Delaunay. First, we compute an unweighted Delaunay triangula- 
tion T Q on M by Algorithm Polyhedral-Delaunay. Then we adjust the weights. The local convexity 
condition is unchanged by an additive constant on all the vertex weights, so T is Delaunay for a vertex 
weighting wq with wq(v) — max u w(u) for all v. Number the vertices v±, . . . , v n in arbitrary order, let vertex 
weighting uii coincide with w on vertices Vj with j < i and with wq elsewhere, and apply algorithm Reweight 
to compute triangulations T\ , Ti , . . . , T„ in turn that are Delaunay for weightings w\ , . . . , w n . Then T = T„ 
is our result. □ 
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Lemma 7.9. Algorithm Polyhedral- Weighted-Delaunay computes a Delaunay triangulation for w in 
0(n 2 log n) time, plus the 0(n 3 S/es) time to compute an unweighted Voronoi diagram. 

Proof. Algorithm Polyhedral-Delaunay costs 0(n 3 S/es) time to compute an unweighted Voronoi dia- 
gram via MMP, plus 0(n) time to convert into an unweighted Delaunay triangulation To. We then apply 
Reweight n times, in time O(nlogn) each, for 0(n 2 logrt) total additional time. 

The output T = T„ is correct by a simple induction using Lemma 17.81 □ 
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